      PROGRAM EXPL2
C***********************************************************************
C                                                                      *
C AUTHOR: YVONNE LUH, GMD, I1.HR                                       *
C                                                                      *
C         POSTFACH 1319                                                *
C         W-5205 ST. AUGUSTIN                                          *
C                                                                      *
C***********************************************************************
C***********************************************************************
C     fortran 90 - Version
C     28/4/92
C.......................................................................
C     EXPLICIT ALGORITHM: STANDARD 2ND ORDER DISCRETIZATION
C     WAVE EQUATION (2 SPACE DIMENSIONS)
C
C     INITIAL BOUNDARY VALUE PROBLEM
C     DDU/DDT=DDU/DDX+DDU/DDY
C     U(X,Y,0)     = F1(X,Y)
C     DU/DT(X,Y,O) = F2(X,Y)
C     U(X,Y,T)     = G(X,Y,T)    (X,Y) AT BOUNDARY OF (0,1)X(0,1)
C.......................................................................
C     YOU CAN TRY DIFFERENT K,H
C     K    STEP IN T-DIRECTION
C     H    STEP IN X- AND Y-DIRECTION
C.......................................................................
C     DIFFERENT EXAMPLES CAN BE CALCULATED
C               COMMON /NEXAMP /NEXAMP
C                1      DSIN(PI*X)*DSIN(PI*Y)
C                            * DCOS(DSQRT(2.0)*PI*T)
C                            * 1000.0D0
C                2      DSIN(20.0*PI*X)*DSIN(20.0*PI*Y)
C                            * DCOS(20.0D0*DSQRT(2.0D0)*PI*T)
C                3      DEXP(X+Y)*DEXP(DSQRT(2.0D0)*T)
C                            * 100000.0D0
C                4      DEXP( 3.0D0*(X+Y))
C                            * DEXP( 3.0D0*DSQRT(2.0D0)*T)
C                5      X**2 + Y**2 + 2.0D0*T**2
C                6      0.0
C     DIFFERENT KINDS OF ERRORS CAN BE CALCULATED
C               COMMON /NRELER  /NRELER
C                1      RELATIVE ERROR    =MAX((SOL - U)/SOL)
C                2      ABSOLUTE ERROR    =MAX (SOL - U)
C***********************************************************************
C***********************************************************************
c
      double precision   u0(:,:), u1(:,:), u2(:,:), tx(:,:), ty(:,:)
      double precision   hu1 (:,:)
      double precision   hu2 (:,:)
      DOUBLE PRECISION   TEND,T,K,H,P
      DOUBLE PRECISION   mxresul, l2resul
      INTEGER            NP,N
      INTEGER            NEXAMP,NRELER
      INTEGER            I, J
c
C.......................................................................
c
C     COMMON BLOCKS
C     NEXAMP             NUMBER OF EXAMPLE
C     NRELER             RELATIVE OR ABSOLUTE ERROR IS CALCULATED
c
      COMMON /NEXAMP /NEXAMP
      COMMON /NRELER /NRELER
c
C***********************************************************************
C***********************************************************************
C     INPUT OF PROBLEM PARAMETERS
C***********************************************************************
C***********************************************************************
C
      READ(5,*) K
      READ(5,*) NP
      READ(5,*) TEND
      READ(5,*) NEXAMP
      READ(5,*) NRELER
C
C***********************************************************************
C***********************************************************************
C     declaration of allocatable arrays
C***********************************************************************
C***********************************************************************
C
      allocate  (u0(0:np,0:np), u1(0:np,0:np), u2(0:np,0:np),
     +           tx(0:np,0:np), ty(0:np,0:np))
      allocate (hu1(0:np,0:np))
      allocate (hu2(0:np,0:np))
C
C***********************************************************************
C***********************************************************************
C     COMPUTATION OF PARAMETERS
C***********************************************************************
C***********************************************************************
C
      N=NP-1
      H=1.0D0/DBLE(NP)
      P=K/H
C
C***********************************************************************
C***********************************************************************
C     calculation of arrays tx and ty
C***********************************************************************
C***********************************************************************
C
!HPF$ INDEPENDENT, LOCAL_ACCESS
      do j = 0, np 
!HPF$ INDEPENDENT, LOCAL_ACCESS
        do i = 0, np
          tx(i,j) = dble (i)
          ty(i,j) = dble (j)
        end do
      end do
        tx = tx * h
        ty = ty * h
c
C***********************************************************************
C***********************************************************************
C     OUTPUT OF PARAMETERS
C***********************************************************************
C***********************************************************************
C
      WRITE(6,*)'*****************************************************'
      WRITE(6,*)'EXPLICIT DISCRETIZATION OF WAVE EQUATION (2D)'
      WRITE(6,*)'*****************************************************'
      WRITE(6,*) ' COMPUTED EXAMPLE: '
      IF      (NEXAMP.EQ.1) THEN
           WRITE(6,*) ' SIN(PI*X)*SIN(PI*Y)*COS(1.414*PI*T) *1000'
      ELSE IF (NEXAMP.EQ.2) THEN
           WRITE(6,*) ' SIN(20*PI*X)*SIN(20*PI*Y)*COS(20*1.414*PI*T)'
      ELSE IF (NEXAMP.EQ.3) THEN
           WRITE(6,*) ' EXP(X+Y)*EXP(1.414*T) * 100000'
      ELSE IF (NEXAMP.EQ.4) THEN
           WRITE(6,*) ' EXP( 3*(X+Y))*EXP( 3*1.414*T)'
      ELSE IF (NEXAMP.EQ.5) THEN
           WRITE(6,*) ' X**2 + Y**2 +2*T**2'
      ELSE IF (NEXAMP.EQ.6) THEN
           WRITE(6,*) ' 0.0 '
      ENDIF
      IF      (Nreler.EQ.1) THEN
           WRITE(6,*) ' relative error is calculated' 
      ELSE IF (Nreler.EQ.2) THEN
           WRITE(6,*) ' absolute error is calculated' 
      ENDIF
      WRITE(6,*) '****************************************************'
      WRITE(6,9000) ' TIME -STEP K = ', K
      WRITE(6,9000) ' SPACE-STEP H = ', H
      WRITE(6,9000) ' RATIO P=K/H  = ', P
C
      IF (P*P.GT.0.5D0) THEN
         WRITE(6,*) ' '
         WRITE(6,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
         WRITE(6,*) ' CFL CONDITION (P*P =< 0.5) NOT FULFILLED'
         WRITE(6,*) ' THIS WILL CAUSE LARGE ERRORS IN COMPUTED SOLUTION'
         WRITE(6,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
         WRITE(6,*) ' '
      ENDIF
C
      WRITE(6,9001) ' NUMBER OF GRID POINTS PER LINE = ',NP
      WRITE(6,9000) ' MAX. T       = ', TEND
 9000 FORMAT(1X,A17,F20.10)
 9001 FORMAT(1X,A35,I5)
      WRITE(6,*) '*****************************************************'
C
C***********************************************************************
C***********************************************************************
C     COMPUTATION OF INITIAL VALUES IN FIRST TWO TIME STEPS
C***********************************************************************
C***********************************************************************
C
      call f1(u0,tx,ty,np)
      call sol(u1,tx,ty,np,k)
      T=K
C
C***********************************************************************
C***********************************************************************
C     TEST OF INITIAL VALUES AT TIME T=K
C***********************************************************************
C***********************************************************************
C
      call difmx(u1,tx,ty,np,k,mxresul)
      call difl2(u1,tx,ty,np,k,l2resul)
      WRITE(6,9020) ' T = ', ' difmx ', ' difl2 '
      WRITE(6,9021) t, mxresul, l2resul
 9020 FORMAT (1X,A8,2x,a10,2x,a10)
 9021 FORMAT (1X,F8.4,2x,D10.4,2x,D10.4)
C
      T=K+K
C
C***********************************************************************
C***********************************************************************
C
   40 CONTINUE
C
C***********************************************************************
C***********************************************************************
C     COMPUTATION OF SOLUTION
C***********************************************************************
C***********************************************************************
C
C     IN INTERIOR
C
C     overlap U1(:,1:1) with HU1
C
C     n = np - 1

      hu1(0:np,1:n) = u1(0:np,2:n+1)
      hu2(0:np,1:n) = u1(0:np,0:n-1)

!HPF$ INDEPENDENT, LOCAL_ACCESS
      do j=1,n
!HPF$ INDEPENDENT, LOCAL_ACCESS
        do i=1,n
         U2(i,j) = P*P* ( U1(i+1,j) + U1(i-1,j) 
     +        + Hu1(i,j) + HU2(i,j)
     +        - 4.0D0 * U1(i,j) )
     +        + 2.0D0 * U1(i,j) - U0(i,j)
       end do
      end do
C
C.......................................................................
C
C     ON BOUNDARY
C
      call g(u2,tx,ty,np,t)
C
C***********************************************************************
C***********************************************************************
C     COMPUTATION AND OUTPUT OF ERROR OF THE COMPUTED SOLUTION
C***********************************************************************
C***********************************************************************
C
      call difmx(u2,tx,ty,np,t,mxresul)
      call difl2(u2,tx,ty,np,t,l2resul)
      WRITE(6,9021) t, mxresul, l2resul
C
C***********************************************************************
C***********************************************************************
C     STORAGE OF NEW SOLUTION
C***********************************************************************
C***********************************************************************
C
      U0 = U1
      U1 = U2
C
C***********************************************************************
C***********************************************************************
C     NEXT TIME STEP
C***********************************************************************
C***********************************************************************
C
      T=T+K
      IF (T.LT.TEND) GOTO 40
C
C***********************************************************************
C***********************************************************************
C     deallocation
C***********************************************************************
C***********************************************************************
c
      deallocate (hu2, hu1)
      deallocate (ty, tx, u2, u1, u0)
c 
C***********************************************************************
C***********************************************************************
C
      END
C
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C
C
      subroutine DIFMX(U,tx,ty,Np,T,mxresul)
C
C     COMPUTES THE MAXIMUM NORM OF THE DIFFERENCE BETWEEN
C     SOL (=SOLUTION OF THE BVP) AND THE VALUES IN U
C
C     DEPENDING ON NRELER  THE RELATIVE ERROR (NRELER=1)
C     OR THE ABSOLUTE ERROR (NRELER=2) IS COMPUTED
C
      integer          np, n
      DOUBLE PRECISION U(0:np,0:np), tx(0:np,0:np), ty(0:np,0:np)
      DOUBLE PRECISION mxresul
      DOUBLE PRECISION uhelp  (0:np,0:np)
      DOUBLE PRECISION uhelp1 (0:np,0:np)
      INTEGER          t
      INTEGER          NRELER
      COMMON /NRELER /NRELER
C
      N = Np-1
C
      call sol(uhelp,tx,ty,np,t)
c      
      uhelp1 = uhelp

      if (nreler.eq.1) then
          where ( uhelp(1:n,1:n) .gt. 1.0d-10 )      
              uhelp1(1:n,1:n) = abs( (uhelp(1:n,1:n)-u(1:n,1:n))
     +                              / uhelp(1:n,1:n) )
          elsewhere
               uhelp1(1:n,1:n) = abs( uhelp(1:n,1:n) - u(1:n,1:n) )
          endwhere
      else if (nreler.eq.2) then
               uhelp1(1:n,1:n) = abs( uhelp(1:n,1:n) - u(1:n,1:n) )
      endif
c
      mxresul = maxval (uhelp1(1:n,1:n))
    
      END
C
C
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C
C
      subroutine DIFl2(U,tx,ty,Np,T,l2resul)
C
C     COMPUTES THE L2 NORM OF THE DIFFERENCE BETWEEN
C     SOL (=SOLUTION OF THE BVP) AND THE VALUES IN U
C
C     DEPENDING ON NRELER  THE RELATIVE ERROR (NRELER=1)
C     OR THE ABSOLUTE ERROR (NRELER=2) IS COMPUTED
C
      integer          np, n
      DOUBLE PRECISION U(0:np,0:np), tx(0:np,0:np), ty(0:np,0:np)
      DOUBLE PRECISION l2resul
      DOUBLE PRECISION uhelp(0:np,0:np)
      DOUBLE PRECISION uhelp1(0:np,0:np)
      INTEGER          t
      INTEGER          NRELER
      COMMON /NRELER /NRELER
C
      N = Np-1
C
      call sol(uhelp,tx,ty,np,t)
c
      uhelp1 = uhelp
      if (nreler.eq.1) then
          where ( uhelp(1:n,1:n) .gt. 1.0d-10 )
              uhelp1(1:n,1:n) = abs( (uhelp(1:n,1:n)-u(1:n,1:n))
     +                              / uhelp(1:n,1:n) )
          elsewhere
               uhelp1(1:n,1:n) = abs( uhelp(1:n,1:n) - u(1:n,1:n) )
          endwhere
      else if (nreler.eq.2) then
               uhelp1(1:n,1:n) = abs( uhelp(1:n,1:n) - u(1:n,1:n) )
      endif
c
      uhelp1(1:n,1:n) = uhelp1(1:n,1:n) ** 2
      l2resul = sum ( uhelp1(1:n,1:n) )
c
      l2resul = dsqrt(l2resul) / dble(n)
c
      END
C
C
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C
C
      subroutine SOL(u,tx,ty,np,T)
C
C**********************************************************************
C
C     CONTINUOUS SOLUTION OF THE DIFFERENTIAL PROBLEM
C
C**********************************************************************
C
      integer          np
      DOUBLE PRECISION u(0:np,0:np), tx(0:np,0:np), ty(0:np,0:np)
      DOUBLE PRECISION t
      DOUBLE PRECISION PI
      INTEGER          NEXAMP
      COMMON /NEXAMP  /NEXAMP
      PI  = DACOS(-1.0D0)
C
C......................................................................
C
      IF (NEXAMP .eq. 1) GOTO 1
      IF (NEXAMP .eq. 2) GOTO 2
      IF (NEXAMP .eq. 3) GOTO 3
      IF (NEXAMP .eq. 4) GOTO 4
      IF (NEXAMP .eq. 5) GOTO 5
      IF (NEXAMP .eq. 6) GOTO 6
C
C......................................................................
C
C     NON OSCILLATING SOLUTION
C
    1 u = dsin(pi*tx)*dsin(pi*ty)*dcos(dsqrt(2.0d0)*pi*t) * 1000.0d0
      GOTO 10
C
C......................................................................
C
C     HIGHLY OSCILLATING SOLUTION
C
    2 u = dsin(20.0d0*pi*tx)*dsin(20.0d0*pi*ty)
     +  * dcos(20.0d0*dsqrt(2.0d0)*pi*t)
      GOTO 10
C
C......................................................................
C
C     SOLUTION INCREASING SLOWLY IN TIME
C
    3 u = dexp(tx+ty) * dexp(dsqrt(2.0d0)*t) * 100000.0d0
      GOTO 10
C
C......................................................................
C
C     SOLUTION INCREASING RAPIDLY IN TIME
C
    4 u = dexp(3.0d0*(tx+ty)) * dexp(3.0d0*dsqrt(2.0d0)*t)
      GOTO 10
C
C......................................................................
C
C     POLYNOMIAL IN X AND Y AS SOLUTION
C
    5 u = tx**2 + ty**2 + 2.0d0*t**2
      GOTO 10
C
C......................................................................
C
C     ZERO SOLUTION
C
    6 u = 0.0D0
      GOTO 10
C
C......................................................................
C
   10 CONTINUE
      END
C
C
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C
C
      subroutine F1(u,tx,ty,np)
C
C**********************************************************************
C
C     INITIAL CONDITION U(X,Y,0) = F1(X,Y)
C
C**********************************************************************
C
      integer          np
      DOUBLE PRECISION u(0:np,0:np), tx(0:np,0:np), ty(0:np,0:np)
      DOUBLE PRECISION PI
      INTEGER          NEXAMP
      COMMON /NEXAMP  /NEXAMP
      PI  = DACOS(-1.0D0)
C
C......................................................................
C
      IF (NEXAMP .eq. 1) GOTO 1
      IF (NEXAMP .eq. 2) GOTO 2
      IF (NEXAMP .eq. 3) GOTO 3
      IF (NEXAMP .eq. 4) GOTO 4
      IF (NEXAMP .eq. 5) GOTO 5
      IF (NEXAMP .eq. 6) GOTO 6
C
C......................................................................
C
C     NON OSCILLATING SOLUTION
C
    1 u = DSIN(PI*tx) * DSIN(PI*ty) * 1000.0D0
      GOTO 10
C
C......................................................................
C
C     HIGHLY OSCILLATING SOLUTION
C
    2 u = DSIN(20.0d0*PI*tx) * DSIN(20.0d0*PI*ty)
      GOTO 10
C
C......................................................................
C
C     SOLUTION INCREASING SLOWLY IN TIME
C
    3 u = Dexp(tx+ty) * 100000.0d0
      GOTO 10
C
C......................................................................
C
C     SOLUTION INCREASING RAPIDLY IN TIME
C
    4 u = Dexp(3.0d0*(tx+ty)) 
      GOTO 10
C
C......................................................................
C
C     POLYNOMIAL IN X AND Y AS SOLUTION
C
    5 u = tx**2 + ty**2
      GOTO 10
C
C......................................................................
C
C     ZERO SOLUTION
C
    6 u = 0.0d0 
      GOTO 10
C
C......................................................................
C
   10 CONTINUE
      END
C
C
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C
C
      subroutine F2(u,tx,ty,np)
C
C**********************************************************************
C
C     INITIAL CONDITION DU/DT(X,Y,0) = F2(X,Y)
C
C**********************************************************************
C
      integer          np
      DOUBLE PRECISION u(0:np,0:np), tx(0:np,0:np), ty(0:np,0:np)
      DOUBLE PRECISION PI
      INTEGER          NEXAMP
      COMMON /NEXAMP  /NEXAMP
      PI  = DACOS(-1.0D0)
C
C......................................................................
C
      IF (NEXAMP .eq. 1) GOTO 1
      IF (NEXAMP .eq. 2) GOTO 2
      IF (NEXAMP .eq. 3) GOTO 3
      IF (NEXAMP .eq. 4) GOTO 4
      IF (NEXAMP .eq. 5) GOTO 5
      IF (NEXAMP .eq. 6) GOTO 6
C
C......................................................................
C
C     NON OSCILLATING SOLUTION
C
    1 u = 0.0d0
      GOTO 10
C
C......................................................................
C
C     HIGHLY OSCILLATING SOLUTION
C
    2 u = 0.0d0
      GOTO 10
C
C......................................................................
C
C     SOLUTION INCREASING SLOWLY IN TIME
C
    3 u = dsqrt(2.0d0) * dexp(tx+ty) * 100000.0d0
      GOTO 10
C
C......................................................................
C
C     SOLUTION INCREASING RAPIDLY IN TIME
C
    4 u = 3.0d0 * dsqrt(2.0d0) * dexp(3.0d0*(tx+ty))
      GOTO 10
C
C......................................................................
C
C     POLYNOMIAL IN X AND Y AS SOLUTION
C
    5 u = 0.0d0
      GOTO 10
C
C......................................................................
C
C     ZERO SOLUTION
C
    6 u = 0.0d0
      GOTO 10
C
C......................................................................
C
  10  CONTINUE
      END
C
C
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C****************************************************************
C
C
      subroutine G(u,tx,ty,np,t)
C
C**********************************************************************
C
C     BOUNDARY CONDITION U(X,Y,T) = G(X,Y,T)
C
C**********************************************************************
C
      integer          np
      DOUBLE PRECISION u(0:np,0:np), tx(0:np,0:np), ty(0:np,0:np)
      DOUBLE PRECISION t
      DOUBLE PRECISION PI
      INTEGER          NEXAMP
      COMMON /NEXAMP  /NEXAMP
      PI  = DACOS(-1.0D0)
C
C......................................................................
C
      IF (NEXAMP .eq. 1) GOTO 1
      IF (NEXAMP .eq. 2) GOTO 2
      IF (NEXAMP .eq. 3) GOTO 3
      IF (NEXAMP .eq. 4) GOTO 4
      IF (NEXAMP .eq. 5) GOTO 5
      IF (NEXAMP .eq. 6) GOTO 6
C
C......................................................................
C
C     NON OSCILLATING SOLUTION
C
    1 u(0,0:np) = dsin(pi*tx(0,0:np)) * dsin(pi*ty(0,0:np))
     +          * dcos(dsqrt(2.0d0)*pi*t) * 1000.0d0
      u(np,0:np) = dsin(pi*tx(np,0:np)) * dsin(pi*ty(np,0:np))
     +          * dcos(dsqrt(2.0d0)*pi*t) * 1000.0d0
      u(1:np-1,0) = dsin(pi*tx(1:np-1,0)) * dsin(pi*ty(1:np-1,0))
     +          * dcos(dsqrt(2.0d0)*pi*t) * 1000.0d0
      u(1:np-1,np) = dsin(pi*tx(1:np-1,np)) * dsin(pi*ty(1:np-1,np))
     +          * dcos(dsqrt(2.0d0)*pi*t) * 1000.0d0
      GOTO 10
C
C......................................................................
C
C     HIGHLY OSCILLATING SOLUTION
C
    2 u(0,0:np) = dsin(20.0d0*pi*tx(0,0:np)) 
     +          * dsin(20.0d0*pi*ty(0,0:np))
     +          * dcos(20.0d0*dsqrt(2.0d0)*pi*t)
      u(np,0:np) = dsin(20.0d0*pi*tx(np,0:np)) 
     +          * dsin(20.0d0*pi*ty(np,0:np))
     +          * dcos(20.0d0*dsqrt(2.0d0)*pi*t)
      u(1:np-1,0) = dsin(20.0d0*pi*tx(1:np-1,0))
     +          * dsin(20.0d0*pi*ty(1:np-1,0))
     +          * dcos(20.0d0*dsqrt(2.0d0)*pi*t)
      u(1:np-1,np) = dsin(20.0d0*pi*tx(1:np-1,np)) 
     +          * dsin(20.0d0*pi*ty(1:np-1,np))
     +          * dcos(20.0d0*dsqrt(2.0d0)*pi*t)
      GOTO 10
C
C......................................................................
C
C     SOLUTION INCREASING SLOWLY IN TIME
C
    3 u(0,0:np) = dexp(tx(0,0:np)+ty(0,0:np)) 
     +          * dexp(dsqrt(2.0d0)*t) * 100000.0d0
      u(np,0:np) = dexp(tx(np,0:np)+ty(np,0:np)) 
     +          * dexp(dsqrt(2.0d0)*t) * 100000.0d0
      u(1:np-1,0) = dexp(tx(1:np-1,0)+ty(1:np-1,0)) 
     +          * dexp(dsqrt(2.0d0)*t) * 100000.0d0
      u(1:np-1,np) = dexp(tx(1:np-1,np)+ty(1:np-1,np)) 
     +          * dexp(dsqrt(2.0d0)*t) * 100000.0d0
      GOTO 10
C
C......................................................................
C
C     SOLUTION INCREASING RAPIDLY IN TIME
C
    4 u(0,0:np) = dexp(3.0d0* (tx(0,0:np)+ty(0,0:np)) ) 
     +          * dexp(3.0d0*dsqrt(2.0d0)*t)
      u(np,0:np) = dexp(3.0d0* (tx(np,0:np)+ty(np,0:np)) ) 
     +          * dexp(3.0d0*dsqrt(2.0d0)*t) 
      u(1:np-1,0) = dexp(3.0d0* (tx(1:np-1,0)+ty(1:np-1,0)) ) 
     +          * dexp(3.0d0*dsqrt(2.0d0)*t) 
      u(1:np-1,np) = dexp(3.0d0* (tx(1:np-1,np)+ty(1:np-1,np)) ) 
     +          * dexp(3.0d0*dsqrt(2.0d0)*t) 
      GOTO 10
C
C......................................................................
C
C     POLYNOMIAL IN X AND Y AS SOLUTION
C
    5 u(0,0:np) = tx(0,0:np)**2 + ty(0,0:np)**2 + 2.0d0*t**2
      u(np,0:np) = tx(np,0:np)**2 + ty(np,0:np)**2 + 2.0d0*t**2
      u(1:np-1,0) = tx(1:np-1,0)**2 + ty(1:np-1,0)**2 + 2.0d0*t**2
      u(1:np-1,np) = tx(1:np-1,np)**2 + ty(1:np-1,np)**2 + 2.0d0*t**2
      GOTO 10
C
C......................................................................
C
C     ZERO SOLUTION
C
    6 u(0,0:np) = 0.0d0
      u(np,0:np) = 0.0d0
      u(1:np-1,0) = 0.0d0
      u(1:np-1,np) = 0.0d0
C
C......................................................................
C
   10 CONTINUE
      END
